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A versatile and numerically inexpensive method is presented allowing the accurate calculation of phase di- 
agrams for bosonic lattice models. By treating clusters within the Gutzwiller theory, a surprisingly good de- 
scription of quantum fluctuations beyond the mean-field theory is achieved approaching quantum Monte-Carlo 
predictions for large clusters. Applying this powerful method to the Bose-Hubbard model, we demonstrate 
that it yields precise results for the superfluid to Mott-insulator transition in square, honeycomb, and cubic lat- 
tices. Due to the exact treatment within a cluster, the method can be effortlessly adapted to more complicated 
Hamiltonians in the fast progressing field of optical lattice experiments. This includes state- and site-dependent 
superlattices, large confined atomic systems and disordered potentials, as well as various types of extended 
Hubbard models. Furthermore, the approach allows an excellent treatment of systems with arbitrary filling fac- 
tors. We discuss the perspectives that allow for the computation of large, spatially-varying lattices, low-lying 
excitations, and time evolution. 
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The recent progress in the realization of optical lattices 
for ultracold atoms offers an ideal testing ground for theo- 
retical models and computational methods. Compared with 
their solid-state counterpart, the foremost advantage of the 
atomic systems is the outstanding controllability of interaction 
strength and lattice parameters. In particular, optical lattices 
allow to study Hubbard models for bosonic particles, where 
an important milestone was the observation of the superfluid 
(SF) to Mott insulator (MI) transition [1,2]. These new exper- 
imental possibilities have stimulated the development of the- 
oretical methods for bosonic systems ranging from density- 
matrix renormalization group (DMRG) techniques and quan- 
tum Monte-Carlo (QMC) methods to mean-field theories. 

In mean-field theories, a single lattice site is decoupled 
from the surrounding lattice, where the fluctuations are de- 
scribed by a mean-field parameter. This strong simplification 
allows nonetheless the qualitative description of strongly cor- 
related phases such as the Mott insulator [2-4]. The decou- 
pling is also achieved in the Gutzwiller method [5-9] where 
the wave function is expanded in local Fock states with in- 
dividual coefficients. Both limits in the phase diagram-the 
Mott state and the superfluid state-can be described (approx- 
imately) by a product of local Gutzwiller states. In fact, it 
can be shown that perturbative mean-field theories and the 
Gutzwiller theory predict equivalent SF-MI transition points 
[6, 10]. The Gutzwiller method is very versatile and can 
be applied to various types of optical lattice systems rang- 
ing from homogeneous lattices to large spatially-varying lat- 
tices systems. It allows both to treat stationary states (e.g. 
Refs. [6, 7, 9, 11]) as well as to perform time evolution (e.g. 
Refs. [9, 12-15]). However, relying on the mean-field decou- 
pling the Gutzwiller theory has similar restrictions as other 
mean-field methods and the phase boundaries can only be cal- 
culated quantitatively. 

Early, it was suggested by Bethe [16], Peierls [17], and 
Weiss [18] that the mean-field method may be extended by 
coupling a cluster of sites (a supercell) with the mean-field 
rather than a single lattice site. Usually, such a cluster is 
formed by a central site and its nearest neighbors. These clus- 
ter mean-field methods were mainly used for the description 
of electrons in solids [19-21] and spin models [22-26] but 



also in the context of dynamical mean-field theories [27]. Re- 
cently, a multi-site mean-field method for the Bose-Hubbard 
model was proposed in Refs. [28-31] based on neglecting 
quadratic terms in the fluctuations. However, due to the 
restriction to relatively small clusters quantum fluctuations 
could not be sufficiently included leading to inaccurate pre- 
dictions of the phase boundary in the SF-MI phase diagram. 
Furthermore, bosonic cluster mean-field theory has been ap- 
plied to dipolar hard-core bosons in Refs. [32, 33]. 

Here, we present a cluster mean-field method for large clus- 
ters based on the well-established Gutzwiller method (Fig. 1). 
As a demonstration, this cluster Gutzwiller method is applied 
to the Bose-Hubbard model in square, honeycomb, and cubic 
lattices. Using periodic boundaries for the cluster, it allows 
for accurate predictions of the SF-MI transition approaching 
quantum Monte-Carlo results for large supercells. We show 
that the approach is suited for the treatment of arbitrary fill- 
ing factors and large-sized lattices. We start by reflecting sev- 
eral aspects of the single-site Gutzwiller method and a general 
description of the supercell method. Subsequently, we apply 
the algorithm to the Bose-Hubbard model in several lattice 
geometries. Finally, details of the method are discussed that 
allow for the computation of large clusters. As an outlook, 




FIG. 1: (a) Within the Gutzwiller theory a single lattice site is cou- 
pled to the mean-field (6.,) of its nearest neighbors (illustrated for a 
square lattice). The lattice sites are expanded in Fock states |n) with 
n = 0, 1, 2, ... particles, (b) In the supercell method, the single site is 
replaced by a cluster of lattice sites and is expanded in the many-site 
Fock basis \N) . 
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we discuss further applications such as the treatment of large 
lattices with non-equivalent sites, the determination of excited 
states or time evolution. 



I. THE SINGLE-SITE GUTZWILLER APPROACH 

Let us start with the standard Gutzwiller method which as- 
sumes that we can write the wave function of the system as 
a product of single-site wave functions. Each site j is repre- 
sented by the Gutzwiller trial wave function \j) = J2 n c n \n) 
[5-9] in the basis of local Fock states \n) with n particles 
(Fig. la). Usually the coefficients c„ are determined by imag- 
inary time evolution [9]. However, as detailed below also a 
self-consistent diagonalization scheme can be used to deter- 
mine the coefficients. Very generally, the wave function of the 
whole system can be written as \j) where \tp) is the wave 
function of all sites except for site j. Common tight-binding 
Hamiltonians can be written as 

H = H^p + Hj + H^pj = H^p + Hj + A^, B", (1) 

where H^ and Hj act only on the respective subsystems. The 
last term represents the coupling between \ip) and |j) and can 
be decomposed in a sum of subsystem operator products. As 
we show in the following, the initial knowledge of ip is not 
required when using a self-consistent loop for the site j. As- 
suming, however, that we know \ip), the Hamiltonian matrix 
of the whole system in the Fock state basis {|n)} is given by 

H mn = (tp\ #V IV') $mn + (H Hj \fl) + 

J2jiP\A%W(m\B?\n), 

where the first term is a constant energy offset. Now, the wave 
function on site j can be obtained by diagonalizing H mn using 
a small Fock basis. For the Bose-Hubbard Hamiltonian 

H w = -jY,b%, + \u^h j {hj-\)- l iY i hj, (3) 

(3,3') 3 3 

with bosonic annihilation (bj) and creation operators on 
site j, on-site interaction U, tunneling matrix element J, and 
chemical potential fi, we obtain for the uncoupled part of the 
Hamiltonian (m\ Hj \n) = \\J n(n— 1) S mn — un 8 mn . The 
coupling H^j between \ip) and \j) is given by 

where the (j 1 ) indicates the summation over all nearest 
neighbors of site j. Using the Fock coefficients c' n of 
the trial wave function on site j', we can obtain (bj,) = 

(ip\ bj, \ip) — J2 n C n c 'n+i^ n + 1- At this point it is clear that 
the Gutzwiller method is equivalent to the mean field treat- 
ment for the Bose-Hubbard model as a single lattice site cou- 
ples only to the average mean field (bj,) and its conjugated. 
In general, the Gutzwiller method has more internal degrees 
of freedom which becomes apparent, e.g., for occupation- 
dependent models [36]. Note that the expression (4) allows 



for the calculation of large-sized inhomogeneous lattices by 
iteratively transversing through the lattice. For a homoge- 
neous lattice with equivalent sites and z nearest neighbors, 
the off-diagonal matrix elements of the coupling are given by 
(n + 1| H^j \n) = —z J sjn + 1 (b) and its conjugated. By 
diagonalizing H mn , a new expectation value of the superfluid 
order parameter (b) = (bj) for the site j can be computed. 
Using a self-consistent loop, the Hamiltonian can be solved 
without initial knowledge of \ip). 

H. THE CLUSTER GUTZWILLER APPROACH 

In the following, the cluster approach is described as an 
extension of the single-site Gutzwiller method. For this it is 
important that the decoupling in equation (2) holds for an ar- 
bitrary subsystem of the lattice rather than a single lattice site. 
Let us therefore replace the single-site state \j) by a super- 
cell cluster |S) with s sites (see Fig. lb). The internal degrees 
of freedom within the supercell allow for quantum fluctua- 
tion that are not covered within the single-site mean-field ap- 
proach. The generalized trial wave function for the cluster is 
given in a many-site Fock basis {|iV)} = {\no, m, 712, ...)} 
reflecting all possible distributions of rii = 0, 1, 2, ... particles 
on the sites i of the supercell. In analogy to the expressions 
(2) and (4), the respective Hamiltonian matrix for the Bose- 
Hubbard model is given by 

H MN = (M\ H S \N) + (V>l K |V) (M\ 4 Q \N) 

= (M\ fl| H - J £ Vj{b)(b) + 6.<St>) \N) , (5) 
jeos 

where the energy offset (ip\H^ \i/j)Smn nas been omitted. 
Here, ijf H is the Bose-Hubbard Hamiltonian (3) where the 
summations are restricted to sites within the supercell S. The 
second term describes the coupling of all sites at the bound- 
ary dS of the cluster with the mean-field (b) . The prefactor v j 
reflects the number of bonds to the mean-field (for Fig. lb 
the corners of the square have Vj = 2 and the other edge 
sites Vj — 1). Note that for more complicated Hamiltoni- 
ans, (ip\A^ \ip) is not reducible to a single mean-field pa- 
rameter (b \) and can contain other internal quantities such as 
(fiibi) obtained from the coefficients Cn- In this case, the 
presented method would differ from the mean-field approach 
in Refs. [28, 29], where fluctuations of quadratic order are ne- 
glected. In equation (5), we further assume a homogeneous 
lattice where all boundary sites couple to the same mean-field 
(b) but it can be easily expanded to super lattices or finite lat- 
tice systems. By solving the eigenvalue problem (5) of the 
cluster (cf. Refs. [28, 29]), we obtain the lowest eigenvector 
| S) =J2n Cn \N), which is used to calculate the mean field 

(b) = (S\ b t \S) = J2 C*mCn (M\ b t \N) . (6) 

M,N 

self-consistently on a target site t. We choose the target site 
t to be the most central site of the supercell. This ensures 
that the mean-field minimally couples to the target site and 
quantum fluctuation can be included as good as possible. Note 
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FIG. 2: (a) The cluster Gutzwiller method applied to the Bose-Hubbard model of a square two-dimensional (2D) lattice. The gray lobe 
corresponds to the single-site Gutzwiller method and the blue lines to supercells with s — 4, 9, and 16 sites (from left to right). Using periodic 
boundary conditions for s — 12 and 16 sites sites (dashed green lines) improve the results significantly. The vertical lines depict the critical 
ratio J/U obtained by mean-field theory (MFT) [3, 4, 6] and quantum Monte-Carlo (QMC). The QMC results (open circles) are taken from 
Ref. [34], The cluster calculation takes / = 7 fluctuations into account (see text), where the red error bar corresponds to / = oo (/ = 8 for 
s > 16) and the black error bar to / = 5. (b) Results for the honeycomb (HC) lattice with three nearest neighbors for s = 4, 12, and 18 
sites (blue) as well as s —18 and 12 sites with periodic boundary conditions (dashed green). The vertical line corresponds to the process chain 
approach (PCA) where the shaded area is the estimated error [35]. (c,d) Finite-size scaling of the critical points for different cluster sizes, 
where circles (diamonds) indicate results without (with) periodic boundary conditions. The scaling parameter A is zero for a single site and 
one for an infinite lattice (see text). For (c) 2D and (d) HC lattices, the predicted critical points for A = 1 match with the QMC and the PCA 
prediction, respectively, within the error bounds. 



that the higher eigenvalues of the Hamiltonian (5) correspond 
to local excitations of the system. 

The feasibility of the cluster Gutzwiller method relies basi- 
cally on three technical aspects, (i) an adequate restriction of 
the many particle basis { | TV) }, (ii) an effective diagonaliza- 
tion procedure, and (iii) an accurate and effective algorithm to 
determine the boundaries in the phase diagram. Before dis- 
cussing the aspects (i)-(iii), we demonstrate the accuracy of 
the cluster Gutzwiller approach by applying it to the Bose- 
Hubbard Hamiltonian (3). 



III. BOSE-HUBBARD PHASE DIAGRAMS 

For two-dimensional square (2D) and honeycomb (HC) lat- 
tices, the resulting phase diagrams for the SF-MI transition 
are shown in Fig. 2. The calculation for s = 1 site cor- 
responds to the Gutzwiller method and therefore is equiv- 
alent to the mean-field (MFT) result with the critical point 
J MFT = (J/J7) crit = 1/(3 + 2y/2)z for the filling factor n = 1 
[3, 4, 6], Here, the number of nearest neighbors is denoted by 
z, where z = 4 for square and z — 3 for honeycomb lat- 
tices. For the square lattice, the prediction J^pj = 0.0429 
using mean-field theory differs substantially from the quan- 
tum Monte-Carlo method [34] and the strong coupling expan- 
sion [38] with Jqmc/sc = 0.0597. By increasing the number 
of sites s in the cluster Gutzwiller method it is possible to 
correct the mean-field result substantially. While for an infi- 
nite number of sites s the mean-field coupling at the boundary 



is negligible and the exact result is obtained, in practice, the 
cluster size is limited to s = 4x4 = 16 sites. Figure 2 demon- 
strates that with an increasing number of sites the results are 
improving. Periodic boundary conditions along one spatial di- 
mensional can improve results for finite-sized lattices signifi- 
cantly, since boundary effects are reduced. Already for a 3 x 3 
cluster in 2D, where the computation is very inexpensive, the 
phase diagram is surprisingly accurate with J|x3 = 0.0559 
at the tip of the Mott lobe. Note that the absolute deviation 
from the numerical exact QMC results (open circles in Fig. 2a) 
is the largest at the tip of the Mott lobe. 

For the accuracy of the cluster method the crucial factor is 
the ratio of internal bonds (within the supercell) to mean-field 
bonds at the boundary. By applying a scaling of the cell size as 
introduced in Ref. [33], the critical J/U for an infinite lattice 
can be interpolated. Here, the scaling with the parameter A = 
Bs/(Bs + Bqs) is performed and plotted Fig. 2c, where Bg 
represents the number of bonds within the cluster and Bqs the 
bonds to the mean field. The scaling parameter also explains 
why periodic boundary conditions (green diamonds) improve 
the numerical results. The predicted value JqEw = 0.0596(4) 
for an infinite lattice (A = 1) matches with the QMC data 
within the standard deviation. 

For the honeycomb lattice (Fig. 2b), the reduced number 
of nearest-neighbors causes an even larger error of the mean- 
field theory value J^pj = 0.0572. The cluster method sub- 
stantially corrects this single-site value and predicts JqEw = 
0.0870(5) for an infinite lattice. This coincides with the pro- 
cess chain approach in Ref. [35] resulting in J P qa = 0.0863 
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FIG. 3: (a) Results for three-dimensional cubic lattices (3D) for su- 
percells with s =4, 12, 18 sites (blue lines) and supercells with 
12 and 18 sites using periodic boundaries (dashed green lines). The 
blue and green line for 18 and 12 sites, respectively, nearly coincide. 
In three dimensions, periodic boundary conditions can be applied in 
two directions for the 3x3x2 cluster (dotted-dashed green line). The 
QMC results (open circles) are taken from Ref. [37]. See Fig. 2 for 
further details, (b) The finite-size scaling for clusters without (cir- 
cles) and with periodic boundary conditions in one (diamonds) and 
two directions (squares). 



within the estimated errors. 

For the three-dimensional cubic lattice (3D), the cluster 
method predicts J c 3 ° w = 0.0350(2) (Fig. 3). The QMC cal- 
culation gives a slightly lower value of J Q 3 ° C = 0.0341 [37]. 
This deviation might be caused by the small edge lengths of 
the clusters in three dimensions. 

The method presented here allows the treatment of systems 
with arbitrary filling factors n. The results for the lowest three 
Mott lobes are shown in Fig. 4. Note that the numerical effort 
does not increase with the filling factors, since only site occu- 
pations rii = n—2, ri+2 are of interest, which is elaborated 
in the next section. 



IV. TREATMENT OF LARGE CLUSTERS 

We now turn back to the description of the numerical algo- 
rithm that allows to treat large cluster sizes. The efficient im- 
plementation allows a numerical inexpensive computation of 
the phase diagram. For a given /i on a todays average (single- 
core) processor, the computation of 9 sites in 2D (5 fluctua- 
tions) takes only few tenths of a second. Here, the method 
also allows to calculate the excitation spectrum. The compu- 
tation time growths drastically with the cluster size, e.g. about 
10 seconds for 16 sites (x 15 for 25 sites and x4 for 6 fluctu- 
ations). 

(i) To apply the method described above, the infinite many- 
particle basis set has to be restricted to a finite but sufficient 
subset { | TV) }. The number of states in {|iV)} grows expo- 
nentially with the particle number and therefore also with 



the number of cluster sites, which limits practically the size 
of the supercell. However, the complexity of the supercell 
problem can be drastically reduced depending on the specific 
Hamiltonian. From the single-site Gutzwiller theory for the 
Bose-Hubbard model, it is known that the MI state with n 
particles is unstable at its boundary only to fluctuations with 
n±l particles; i.e., only zero, one and two particle Fock 
states have to be taken into account for the MI with a filling 
of n — 1 [6], This is, however, not completely true for the 
supercell method which has internal degrees of freedom. It 
turns out that local fluctuations with rii = ±2 have an effect, 
if rather small, whereas higher local particle number fluctua- 
tions are extremely small and can also be neglected (see er- 
rorbars in Figs. 2 and 3). However, we can indeed restrict 
the total particle number to N and N ± 1. Moreover, at the 
phase boundary also states with all sites simultaneously fluc- 
tuating are unlikely to be occupied. In practice, more than 
/ = J2i I n i ~ n I < 5 fluctuations for smaller clusters and 
f < 7 for larger clusters (s > 16) hardly change the results 
of the phase diagram (see error bars in Fig. 2 and Fig. 3). In 
the vicinity of the boundary the actual value of the superfluid 
order parameter is influenced by this constraints but not the 
criticality. 

(ii) For the diagonalization of the coupled supercell prob- 
lem Lanczos-based algorithms for sparse matrices can be used 
to obtain the lowest eigenvector [39]. In fact, the structure of 
the Hamilton matrix has to be created only once for a given 
number of lattice sites. As usually consecutive points of a 
phase diagram are calculated, excellent guesses for both ip 
and the eigenvector can be provided speeding up the Lanc- 
zos diagonalization immensely. Furthermore, symmetries of 
the cluster can be used to build a symmetrized basis set, e.g., 
for a 4 x 4 cluster with periodic boundaries the C-2 and the C4 
symmetry allow a reduction of the basis length by a factor of 
8. While here we restrict ourselves to the ground state, this 
optimized method is also suited to compute a large number of 
excited states. 
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FIG. 4: Mott insulator phase boundaries for filling factors n = 1-3 
obtained by the supercell approach, which allows for the calculation 
of arbitrary filling factors. The results for the cubic (3D), the square 
(2D), and the honeycomb lattice (HC) are plotted for the largest cell 
sizes with periodic boundary condition along one direction as a func- 
tion of zJ/U. In this unit the mean-field boundary (gray) is indepen- 
dent of the lattice geometry. 
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(iii) A straightforward implementation would calculate all 
points in the fi/U - J/U plane of the SF-MI phase diagram. 
In addition, the self-consistent loop convergences only slowly 
directly at boundary due to the sudden change of the super- 
fluid order parameter ip = (b), which would require many 
iterations to determine the phase boundary accurately. How- 
ever, it is much more convenient to rely only on the fact that 
the algorithm converges monotonically. Starting with an ini- 
tial guess (f t h (threshold value), a single iteration allows to de- 
termine whether the exact value of ip is smaller or greater than 
tpth. As the value of ip has a jump at the boundary, a small but 
finite (^ t h, such as 10 -6 , allows therefore to determine whether 
a given point is in the SF or the MI phase. Applying a binary 
search algorithm with only 20 iterations for a given value of ^ 
allows to determine the critical value with a relative precision 
of2- 20 «l(r 6 . 



V. CONCLUSIONS 

We have presented an intuitive cluster mean-field method 
based on the Gutzwiller theory where the single-site Fock 
state is replaced by a supercell. Using large clusters, we 
demonstrate for the Bose-Hubbard Hamiltonian that this 



method allows accurate results for various lattice geometries 
and arbitrary filling factors, in particular, for 2D cubic and 
honeycomb lattices. The approach can be adapted to various 
types of Hamiltonians and is numerical inexpensive. An in- 
trinsic advantage of the method is that it can be easily ap- 
plied to very large lattices with inequivalent lattice sites [29], 
such as disorder potentials or macroscopic confined systems. 
To achieve this, the supercell centered at site x can be iter- 
atively moved through the lattice where each time the mean 
field ip x = (b x ) at the target site x is updated. Another ad- 
vantage is that the method can be used to compute the local 
excitation spectrum. Furthermore, it also allows to perform 
correlated time evolution. For each short time step t i+1 and 
for each lattice site x t , the exact time-evolution can be per- 
formed within the supercell centered at site x t using the time- 
dependent mean-field boundary ip x (U). This determines the 
expectation values ¥> x (*i+i) for the next time step. 
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